/*
  This file is part of code_saturne, a general-purpose CFD tool.

  Copyright (C) 1998-2025 EDF S.A.

  This program is free software; you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation; either version 2 of the License, or (at your option) any later
  version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  details.

  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
  Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

/*----------------------------------------------------------------------------*/

#include "base/cs_defs.h"

/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*----------------------------------------------------------------------------
 * Local headers
 *----------------------------------------------------------------------------*/

#include "bft/bft_error.h"
#include "bft/bft_printf.h"

#include "base/cs_assert.h"
#include "base/cs_base.h"
#include "base/cs_field.h"
#include "base/cs_field_pointer.h"
#include "base/cs_log.h"
#include "base/cs_math.h"
#include "base/cs_mem.h"
#include "base/cs_time_step.h"
#include "mesh/cs_mesh.h"
#include "mesh/cs_mesh_location.h"
#include "mesh/cs_mesh_quantities.h"
#include "base/cs_numbering.h"
#include "base/cs_parall.h"
#include "base/cs_physical_constants.h"
#include "pprt/cs_physical_model.h"

#include "atmo/cs_atmo.h"
#include "atmo/cs_air_props.h"
#include "atmo/cs_atmo_aerosol_ssh.h"
#include "atmo/cs_atmo_chemistry.h"
#include "atmo/cs_intprf.h"

/*----------------------------------------------------------------------------
 *  Header for the current file
 *----------------------------------------------------------------------------*/

#include "atmo/cs_atmo_kinetic_rates.h"

/*----------------------------------------------------------------------------*/

BEGIN_C_DECLS

/*============================================================================
 * Prototypes for external Polyphemus functions generated by SPACK.
 *============================================================================*/

void
cs_ext_polyphemus_kinetic_1(int       nr,
                            double    rk[],
                            double    temp,
                            double    xlw,
                            double    press,
                            double    azi,
                            double    att,
                            int       option_photolysis);

void
cs_ext_polyphemus_kinetic_2(int       nr,
                            double    rk[],
                            double    temp,
                            double    xlw,
                            double    press,
                            double    azi,
                            double    att,
                            int       option_photolysis);

void
cs_ext_polyphemus_kinetic_3(int       nr,
                            double    rk[],
                            double    temp,
                            double    xlw,
                            double    press,
                            double    azi,
                            double    att,
                            int       option_photolysis);

void
cs_ext_polyphemus_kinetic_4(int       nr,
                            double    rk[],
                            double    temp,
                            double    xlw,
                            double    press,
                            double    azi,
                            double    att,
                            int       option_photolysis);

END_C_DECLS

/*=============================================================================
 * Additional doxygen documentation
 *============================================================================*/

/*!
  \file cs_atmo_kinrates.cpp
        Computation of reaction rates for atmospheric chemistry
 */

/*! \cond DOXYGEN_SHOULD_SKIP_THIS */

/*=============================================================================
 * Local macro definitions
 *============================================================================*/

/*============================================================================
 * Type definitions
 *============================================================================*/

typedef void
(cs_atmo_kinrates_func_t) (int       nr,
                           double    rk[],
                           double    temp,
                           double    xlw,
                           double    press,
                           double    azi,
                           double    att,
                           int       option_photolysis);

/*============================================================================
 * Global variables
 *============================================================================*/

/*============================================================================
 * Private function definitions
 *============================================================================*/

/*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */

BEGIN_C_DECLS

/*=============================================================================
 * Public function definitions
 *============================================================================*/

/*----------------------------------------------------------------------------*/
/*!
 * \brief  Compute reaction rates for atmospheric chemistry.
 */
/*----------------------------------------------------------------------------*/

void
cs_atmo_kinetic_rates_compute(void)
{
  // Initialization

  const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
  const cs_real_3_t *cell_cen = cs_glob_mesh_quantities->cell_cen;

  const cs_atmo_option_t *at_opt = cs_glob_atmo_option;
  const int meteo_profile = at_opt->meteo_profile;

  cs_fluid_properties_t *fluid_props = cs_get_glob_fluid_properties();
  cs_real_t rair = fluid_props->r_pg_cnst;

  cs_real_t *cvar_totwt = nullptr, *cpro_liqwt = nullptr;
  cs_real_t *cpro_met_temp = nullptr, *cpro_met_p = nullptr;
  cs_real_t *cpro_tempc = nullptr, *cpro_met_qv = nullptr;
  cs_real_t *crom = nullptr;

  if (meteo_profile == 2) {
    cpro_met_temp = cs_field_by_name("meteo_temperature")->val;
    cpro_met_p = cs_field_by_name("meteo_pressure")->val;
    cs_field_t *f_met_qv = cs_field_by_name_try("meteo_humidity");
    if (f_met_qv != nullptr)
      cpro_met_qv = f_met_qv->val;
  }

  cs_atmo_chemistry_t *at_chem = cs_glob_atmo_chemistry;
  const int iphotolysis = (at_chem->chemistry_with_photolysis) ? 1 : 2;

  int atmo_model_flag = cs_glob_physical_model_flag[CS_ATMOSPHERIC];
  if (atmo_model_flag >= 1) {
    crom = CS_F_(rho)->val;
    cpro_tempc = cs_field_by_name("real_temperature")->val;
  }

  if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] >= 2) {
    cvar_totwt = cs_field_by_name("ym_water")->val;
    cpro_liqwt = cs_field_by_name("liquid_water")->val;
  }

  // Computation of zenith angle
  cs_real_t qureal = (cs_real_t)at_opt->squant;             // julian day
  cs_real_t ut_time =   (cs_real_t)at_opt->shour            // UTC time
                      + (cs_real_t)at_opt->smin / 60.0
                      + (cs_real_t)at_opt->ssec / 3600.0;

  int idtvar = cs_glob_time_step_options->idtvar;
  if (idtvar == 0 || idtvar == 1)
    ut_time += cs_glob_time_step->t_cur / 3600.0;

  cs_real_t albedo = 0;    // albedo, unused here
  cs_real_t za = 0;        // zenith angle (rad)
  cs_real_t dlmuzero = 0;  // cos of zenith angle
  cs_real_t omega = 0;
  cs_real_t fo = 0;        // solar constant, unused here

  cs_atmo_compute_solar_angles(at_opt->latitude,
                               at_opt->longitude,
                               qureal,
                               ut_time,
                               0, /* no Sea */
                               &albedo,
                               &za,
                               &dlmuzero,
                               &omega,
                               &fo);

  cs_real_t azi = cs::abs(za*18.0/cs_math_pi);  // zenith angle (degree)

  /* Computation of kinetic rates in each cell
     ----------------------------------------- */

  // Note: Photolysis should be cut in SPACK and not here
  //       even if the azimuthal angle is > 90

  cs_real_t t_cur = cs_glob_time_step->t_cur;

  int met_1d_nlevels_t = at_opt->met_1d_nlevels_t;
  int met_1d_ntimes = at_opt->met_1d_ntimes;
  const cs_real_t *z_temp_met = at_opt->z_temp_met;
  const cs_real_t *time_met = at_opt->time_met;
  const cs_real_t *hyd_p_met = at_opt->hyd_p_met;
  const cs_real_t *qw_met = at_opt->qw_met;

  cs_atmo_kinrates_func_t  *kinrates_func;

  switch(at_chem->model) {
  case 1:
    kinrates_func = cs_ext_polyphemus_kinetic_1;
    break;
  case 2:
    kinrates_func = cs_ext_polyphemus_kinetic_2;
    break;
  case 3:
    kinrates_func = cs_ext_polyphemus_kinetic_3;
    break;
  case 4:
    kinrates_func = cs_ext_polyphemus_kinetic_4;
    break;
  default:
    kinrates_func = nullptr;
  }

  cs_assert(kinrates_func!= nullptr);

  const cs_lnum_t nrg = at_chem->n_reactions;
  cs_real_t *reacnum = at_chem->reacnum;

  // Loop on cells

  int n_loc_threads = cs_parall_n_threads(n_cells, CS_THR_MIN);

  double *rk_buf;
  size_t rk_buf_t_size = cs_numbering_simd_size(nrg, sizeof(cs_real_t));

  CS_MALLOC(rk_buf, rk_buf_t_size*n_loc_threads, double);

  #pragma omp parallel num_threads(n_loc_threads)
  {
#if defined(HAVE_OPENMP)
    const int t_id = omp_get_thread_num();
#else
    const int t_id = 0;
#endif

    cs_lnum_t t_s_id, t_e_id;
    cs_parall_thread_range(n_cells, sizeof(double), &t_s_id, &t_e_id);
    double *rk = rk_buf + t_id*rk_buf_t_size;

    for (cs_lnum_t c_id = t_s_id; c_id < t_e_id; c_id++) {

      // Temperature and density

      cs_real_t temp, press;

      // Dry or humid atmosphere
      if (atmo_model_flag >= 1) {
        temp = cpro_tempc[c_id];
        press = crom[c_id] * rair * temp;   // ideal gas law
      }

      // Constant density with a meteo file
      else if (meteo_profile == 1) {

        cs_real_t zent = cell_cen[c_id][2];  // Z coordinate of the cell

        // Hydrostatic pressure
        press = cs_intprf(met_1d_nlevels_t,
                          met_1d_ntimes,
                          z_temp_met,
                          time_met,
                          hyd_p_met,
                          zent,
                          t_cur);

        // Temperature
        temp = cs_intprf(met_1d_nlevels_t,
                         met_1d_ntimes,
                         z_temp_met,
                         time_met,
                         z_temp_met,
                         zent,
                         t_cur);

      }

      else {
        press = cpro_met_p[c_id];
        temp = cpro_met_temp[c_id];
      }

      temp += cs_physical_constants_celsius_to_kelvin;

      // Specific humidity

      cs_real_t hspec;

      // Humid atmosphere
      if (atmo_model_flag >= 2)
        hspec =   (cvar_totwt[c_id]-cpro_liqwt[c_id])
                / (1.0-cpro_liqwt[c_id]);

      // Constant density or dry atmosphere with a meteo file
      else if (meteo_profile == 1)
        hspec = cs_intprf(met_1d_nlevels_t,
                          met_1d_ntimes,
                          z_temp_met,
                          time_met,
                          qw_met,
                          cell_cen[c_id][2],
                          t_cur);

      else if (cpro_met_qv != nullptr)
        hspec = cpro_met_qv[c_id];

      // Dry
      else
        hspec = 0;

      // Compute kinetic rates
      kinrates_func(nrg, rk, temp, hspec, press, azi, 1., iphotolysis);

      // Store kinetic rates
      for (cs_lnum_t ii = 0; ii < nrg; ii++) {
        reacnum[ii*n_cells + c_id] = rk[ii];
      }

    } // End of loop on cells

  } /// End of OpenMP section

  CS_FREE(rk_buf);
}

/*----------------------------------------------------------------------------*/

END_C_DECLS

